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Abstract 

While enjoying demonstrated improvement in accuracy, efficiency, and robustness over 
existing schemes, the Advection Upstream Splitting Scheme (AUSM)[1,2] has been found 
to have some deficiencies in extreme cases. This paper describes recent progress towards 
improving the AUSM while retaining its advantageous features. The new scheme, termed 
AUSM + , features: (1) Unification of velocity and Mach number splittings reported in [1,2], 
(2) Exact capture of a single stationary shock, and (3) Improvement in accuracy. 

In this paper, we lay out a general construction of the AUSM + scheme and then focus 
on the analysis of the scheme and its mathematical properties, heretofore unreported. We 
prove monotonicity and positivity, and give a CFL-like condition, for first and second 
order schemes and for generalized curvilinear co-ordinates. Finally, results of numerical 
tests on many problems axe given to confirm the capability and improvements on a variety 
of problems including those failed by prominent schemes. 


1. INTRODUCTION 

Seeking an accurate numerical scheme for capturing shock and contact discontinuities, 
with minimal n um erical dissipation and oscillations, has been a lasting challenge to the 
computational fluid dynamicist as well as numerical analyst since the advent of comput- 
ers. The 1980s witnessed an explosive interest and research in upwind schemes for their 
capability of achieving high accuracy over a wide range of problems described by Euler 
or Navier-Stokes equations. Today, upwind schemes undoubtedly have become the main 
spatial discretization techniques adopted in nearly all major codes. 

While much has been said in the literature about the successes and fails of various 
schemes (see for example a recent survey by Roe [3]), Quirk[4] recently added an interest- 
ing catalogue of situations that deny several current Riemann solvers success, serving to 


increase the general level of awareness to their limitations. This report will, however, only 
comment on some of these methods to the extent that is relevant to the motivations and 
issues concerning the the development of a new flux scheme. 

It is well known that upwind schemes are generally classified in two categories, either 
flux-vector or flux-difference splittings. Flux-vector splittings (FVS) have proved to be 
a simple and useful technique for arriving at upwind differencing and are pre-eminently 
suited for use in implicit schemes. The Steger- Warming splitting [5] utilizes the com- 
plete wave struture (the eigenvectors of the Jacobian matrix — three waves in ID Euler 
equations), and the homogeneous property at a given state. This splitting, achieved by 
grouping positive and negative eigenvalues, is nondifferentiable and leads to glitches when 
eigenvalues change sign. Fixes have been proposed to remedy this deficiency, but at the 
expense of introducing excessive numerical diffusion. 

The Van Leer splitting [6] uses only two waves to achieve splitting. The splitting is 
made to be differentiable and gives accurate solution for many steady and unsteady inviscid 
problems. Our recent experiences, to be reported separately, further reaffirm that the VLS 
is quite an admirable solver for inviscid problems, in particular those involving strong 
shocks and their interactions. Unfortunately the intent for making a smoother transition 
at the sonic points brings an unacceptable amount of diffusion at low speed, attaining 
maximum error at M = 0. As a result, the boundary layer which has nearly vanishing 
transverse velocity is significantly broadened, yielding incorrect pressure and temperature 
distributions [1,7]. Owing to the same cause, significant error also appears for stationary 
or slowly moving contact discontinuity. (See Fig. 5.) A strategy for eliminating the above 
diffusion has been proposed in [8] in which each mass flux component is made to vanish as 
\M\ = 0. This, however, turns out to be too strong a requirement for designing numerical 
fluxes — a serious stability problem crops up. 

Flux difference splittings (FDS) are generally acknowledged to be accurate for both 
Euler and Navier-Stokes calculations, although they are substantially more expensive than 
the FVS to compute. The Roe splitting [9] perhaps is most popular of the FDS; nearly 
every major code offers an option to use the Roe splitting in one form or another. A major 
complaint about the Roe scheme is its need to inject an entropy fix to eliminate certain 
nonphysical solutions. Unfortunately this fix should be tuned between problems, inevitably 
introducing superfluous diffusivity to the viscous solution in exchange for stability. In 
addition to the failings catalogued by Quirk[4], the Roe splitting is found to be unstable 
for unsteady rarefaction-rarefaction waves. 

The Osher-Soloman splitting [10], save its entropy-satisfying property, is found not 
totally immune from difficulties. Its solution may be path-dependent on the natural or 
reverse order chosen. One example where the Osher-Solomon scheme fails is again the 
bl un t body problem where the natural order is found to yield an anomalous solution, 
although different from Roe’s solution. 

A first glance of the above summary may suggest a bleak conflict — the simplicity 
of the FVS comes at a price of reduced accuracy due to numerical diffusion, while the 
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FDS attains increased accuracy at a cost of decreased efficiency. The current research 
is motivated by the desire to combine the efficiency of FVS and the accuracy of FDS to 
arrive at ways that do away with the problem of numerical diffusion with only a small 
(if any) increase in complexity. To be practically useful, such a scheme should hold ro- 
bustness/stability, in addition to accuracy and efficiency, for a wide range of problems 
Euler and Navier-Stokes equations, ideal and nonequilibrium gass, and steady and un- 
steady flows. Our recent attempts towards deriving a scheme meeting these goals have 
proven to be quite fruitful, resulting in so-called the AUSM [1,2] and its derivatives such 
as AUSMDV [11] and the present AUSM+, and the HUS [12,13]. Both AUSMs and HUS 
are based upon rather different design concepts and the latter, reported earlier in [12,13] by 
Coquel and Liou, will be further elaborated with a much detailed theoretical development 
and construction in a forthcoming paper [14]. Another product deserving serious attention 
as well is called AUSMDV and is given with full details in [11]. The present paper is a 
complement to the above research, but addresses different topics of interest. 

Several efforts have been attempted to improve the original Van Leer scheme, in 
particular by Hanel et al [15,16]. Van Leer, recently based on their approach [16], has made 
a significant improvement [17], in which the temperature distribution of a hypersonic conic 
flow is predicted accurately. Unfortunately, pressure glitch is observed in the calculation. 
Thus, the effort for searching improved numerical flux functions is still continuing, striving 
for the goal of arriving towards the so-called “ultimate” scheme. 

Unlike the conventional flux schemes such as flux- vector and flux-difference splittings, 
we begin the development of the AUSM schemes by separating the convective and pressure 
fluxes. Appropriate polynomial functions for the nonlinear waves ((u ± a) in ID case) are 
used to represent their interactions, leading to representations of the convective velocity 
and the pressure at the cell interface. Finally the scalar quantity (p, pu, pH ) at the upwind 
location is properly advected by this interface velocity, hence the scheme is coined as the 
Advection Upstream Splitting Method, AUSM for short. The new flux scheme has been 
shown to yield improved accuracy, efficiency, and robustness by us and independently by 
others. 

Most recently, the notion of this splitting procedure has been adopted by several au- 
thors [18,19] for various motivations. Jameson [18] rewrote the splitting in the form of 
artificial dissipation for the central-differencing setting and chose a different set of poly- 
nomial representations for splitting M and p. Unfortunately, an adjustable and problem- 
dependent constant has been inserted in the formula, thereby reverting back to complaints 
that upwind differencing is meant to avoid. For usual problems, a nonvanishing constant 
has to be set, leading to excessive smearing in both contact and shock discontinuities. 
Halt and Agarwal [19] have experimented with a variant in which the pressure work term 
is separated out from the total enthalpy, as first proposed by Steger [5,20]. The immediate 
drawback of doing this is that the total enthalpy is no longer guaranteed to be conserved 
when the situation so demands (see cases and related results in [21]). The conservation 
of total enthalpy, not only having theoretical justification, is also proved to be critical 
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in viscous calculations. Indeed, this is the very idea of Hanel [15] for retaining the total 
enthalpy in Van Leer’s FVS. 

Even as the original scheme AUSM enjoys remarkable success, its drawbacks also 
surface. In a continuing search for a near-perfect (if it exists) numerical flux scheme, we 
report in this paper a recent progress towards this end. We aim at deriving an improved 
flux formula by addressing some fundamental properties by way of rigorous mathematical 
derivation. Consequently, this new-revised scheme, termed AUSM + , can be shown either 
analytically or numerically, to satisfy the following properties: (1) Exact resolution of a 
single stationary normal shock or contact discontinuity, (2) Positivity condition, and (3) 
Entropy condition. Property (1) is important not only in shock resolution of inviscid flows, 
but also in boundary /shear- layer resolution of viscous flows. Property (2) proves to be 
closely related to the robustness of the scheme for high speed flows and species calculation. 
And property (3) selects the physically correct solution by rejecting the expansion shocks. 
While these are desirable properties, none of the existing schemes are known to possess 
all of them. For example, properties (2) and (3) are violated by the Roe splitting, and (1) 
and (2) by the Osher splitting. In practice, some “fixes”, of little mathematical or physical 
justification, have been inserted to prevent the scheme from failing. In the case of the Roe 
scheme, a heavy dose of entropy fix can be prescribed to cure the “carbuncle” phenomena, 
which is believed to be caused primarily by its failure to satisfy property (3). 

Like its predecessor, the AUSM + is simple conceptually with close connection to 
physical interpretation; it is also simple algebraically, leading to greater computational 
efficiency, for it requires only 0(n) operations for evaluating fluxes, n being the dimension 
of the flux. 

In this paper, we will give the details of the development and analysis of the AUSM + 
scheme for ideal-gas flows. In a separate paper [21], we will include non-equilibrium flows. 
This paper is organized as follows. In section 2 we give a detalied derivation of the proposed 
scheme and prove some relevant properties. In section 3 we analyze the numerical flux 
developed in the previous section and compare it with other flux schemes. Also we present 
proof of monotonicity and positivity, and the associated CFL-like condition. In section 4 
we give formulas for extension to higher-order accuracy and present similar proofs as in 
section 3. In section 5 an extention to the generalized curvilinear coordinates is given. 
In section 6 we discuss the boundary conditions used in the calculations. Test cases are 
shown in section 7 together with discussion and comparisons with other schemes. Finally 
we conclude the paper with a brief remark. 

2. FORMULATION 


Consider as an initial- value problem the one- dimensional system of conservation laws 
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for ideal-gas flows: 


_ +g -=0 , <>0 , -oc<x<oo, (!) 

U(x,0) = Uo(x), 

where U = ( p , pu, pet) T belongs to a phase space U € 'RZ, the inviscid flux F = (pu, pu 2 + 
p,puh t ) T denotes a smooth mapping F : U -* K p , and the specific total energy e t = 
e + u 2 / 2 = ht-p/p. 

For the numerical solution of (1), we shall consider piecewise constant approximations 
U” +1 defined by the explicit 3-point scheme in conservation form 

U" +1 = UJ - A(F” +1/2 - F"_ 1/2 ), n € N, j € Z, (2) 

where A = At /Ax, At and Ax being respectively the time and space steps. Here, the 
numerical flux defined by 

*? h/ ,-F <+1 / 1 (U?,UJ +1 ). 5 (3) 

is assumed to be a Lipschitz continuous function and satisfy the consistency condition: 

F i+1/2 (U,-,U i ) = F(U J ). (4) 


In the finite-volume formulation (2)-(3), the differences among all numerical schemes 
lie essentially in the definition of the numerical flux F y+ 1/2 evaluated at the cell interface. 
It is obvious that different numerical fluxes will satisfy different sets of criteria such as 
positivity, entropy condition, differentiability, etc. and are expected to yield various degrees 
of satisfaction in performance. 

As a first step in the formulation of the AUSM family of schemes, we recognize the 
convection and acoustic waves as two physically distinct processes, and write the inviscid 
flux as a sum of the convective and pressure terms: 

F = F (c) + P, (5a) 


where 


F (c) = c$ = I 


u j pu j , if c — u 

\phtj 

( P a \ 

M I pua I , if c = M, 

\ph t aj 


(56) 


§ In this paper we always use the half-integer subscripts “j + 1/2”, j € Z, to denote 
“numerical” function f : U x IA — t W 6 1Z? . It should not be confused with a variable 
having an integer subscript “ 7 ”, which indicates a function / : U — *• U € 'RZ . 
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and 



Here the convective flux contains the passive scalar quantities in $ convected by the 
velocity c. Depending on whether the convective speed c is chosen to be either u or M, the 
vector $ has the corresponding contents. The pressure flux P contains solely the pressure 
term. 

One immediate question to be posed is which of the two convective speed formulas, 
respectively referred to as the U- and M-splittings, is preferred. In [2], we stated that the 
U-splitting was found to be more robust than the M-splitting for the shock-tube problem 
since a larger time step is allowed at start, while solutions for several other test problems 
gave no significant differences. On the other hand, the former smears the stationary 
contact discontinuity since the interface velocity does not vanish, while the latter holds 
the contact discontinuity exactly. This ambiguity actually suggests that a unified formula 
should be pursued, which indeed leads to a major discovery in our continuing effort towards 
improving the numerical flux. The net result is a unified formula that in addition leads to 
an exact capturing of a single stationary shock. Moreover, due to a judicious choice of a 
common speed of sound for defining the interface Mach number, the new numerical flux 
can give rise to an exact capture of a single stationary shock. 

As a preliminary, we shall express the numerical convective flux F^ 1 y 2 at the interface 
j + 1/2 straddling the j-th and (j + l)-th cells as 

F 5+i/2 = M j+l/2®j+l/2- (6) 

We stress now that, as will be seen later, it is immaterial in our development whether 
we write it in terms of u or M for the unified formula. As the M-splitting may be more 
familiar, we will adopt it in this section and change to the U-splitting in the next section 
for a cleaner algebraic manipulation. 

To define the numerical fluxes, two steps are taken. First, we define the interface 
velocity and pressure by considering the nonlinear fields. Mathematically, we propose to 
separately deal with the genuinely nonlinear fields (u + a, u — a) and the linearly degenerate 
field ( u ). In this approach, the coupling between velocity and pressure, as described by 
the characteristic relationships, is honored via the genuinely nonlinear fields, which will be 
described in detail in §2.2. In the sequel we respectively give detailed formulas for the major 
elements involved in F j+x/ 2 : (1) definition of $/+ 1 / 2 ? (2) definition of (Mj+x/2-> Pj+i/ 2 )i 
and (3) definition of a J+1 /2- 

2.1 Definition of 4>y + x /2 
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The scalar quantities in $y-|- 1/2 are convected by the linear field via the interface 
velocity My + i/ 2 using a simple upwinding: 


* i + 1/2 = 


*i+i» 


if Mj+i/2 > 0, 
otherwise. 


(7) 


2.2 Definition of (My+i/ 2 ,p 7 +i/ 2 ) 

Anticipating contributions from a j” and u j + 1” states to the interface Mach number 
M j+1 / 2 , let us write the mapping m :U xU —*U € 

Mj+i /2 = m(Mj , M j+ 1 ). (8a) 

Specifically we write My+j/ 2 as a sum of two individual components, 

M j+1/2 = M/ + M- +1 , (86) 

where the superscripts and u — ” are understood to be associated with the right- and 
left-running waves. 

For ID conservation laws, the nonlinear characteristic equations are: 

dp ± padu = 0, along — = u ± a. (9) 

This simply suggests that the velocity and pressure are closely coupled and both quantities 
may be thought of propagating locally at the speeds tt±a. For \M\ < 1, the two waves move 
in the opposite directions and interact with each other. Hence the interface velocity and 
pressure, composed of the interaction of these two waves, are constructed using (uia) as 
basis functions for polynomial expansions. This suggests the following expansion in terms 
of (M ± 1): 

M ± = \{M ± 1), for \M\ < 1. (10) 

This in fact is a form of eigenvalue expansion, just utilizing the genuinely nonlinear fields. 
Then from (8b), the interface Mach number is defined, as \M\ < 1, by 

Mj+1/2 = i[(Af + 1), + (M - l)y+i). (11) 

Note that as in Van Leer’s splitting, the present method heavily makes explicit use of the 
eigenvalues of the nolinear wave, Mil. Using the eigenvalues as a basis for expressing 
the numerical fluxes is quite common in the upwind formulation. For example, it comes 
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naturally in the Steger- Wanning and Roe fluxes as the Jacobian matrix is explicitly ex- 
pressible in terms of its eigenvalues. Specifically, Steger- Warming splitting [5] begins by 
grouping F as: 


F = 


P_ 

2 7 


(u + a) 


( 1 

I u + a 

\H + ua 


+ (u 



-I- 2(7 - l)u 



( 12 ) 


To remove non-differentiability of (10) when sign changes, Van Leer [6] chose differ- 
entiable, second-order polynomials: 

M ± = ±J(M±1) 2 . (13) 

It is evident that the above split formula (13) results in an unsymmetric distribution of 
signals propagated by the characteristic speeds ( M i 1) as M ^ 0 and becomes the simple 
(symmetric) average as M tends to zero, leading to the formula used commonly in the 
pressure-based codes for incompressible flows. 

In the present study, we present a more general footing for designing the split formulas. 
First we introduce the following: 

Proposition 2.1: Let the respective split Mach numbers M ± be chosen such that they 
hold the following properties: 

[Ml]: M + + M~ = M, for consistency. 

[M2]: M + > 0 and M~ < 0. 

[M3]: M ± are monotone increasing functions of M. 

[M4]: M + (M) = i.e., a symmetry property. 

[M5]: M + = M as M > 1; M~ = M as M < —1. 

[M6]: M ± are continuously differentiable. 

Consequently we assert the following necessary results. 

Lemma 2.1: Mj + ij 2 € [Mj,Mj+i], j € Z. 

Proof. Let D = (M j+ 1/2 - Mj)(Mj +1 f 2 ~ M j+ 1 ), substituting (86) and [Ml] in D gives 

D = -(M- - - M* +J ) < 0 

where the product of two parentheses is nonnegative by virtue of [M3]. 

Remark: Owing to [M2], [M3], and [M5], it is easy to see 0 < M + < 1 and 0 > M~ > -1 
as \M\ < 1. 

Definition of M: The split Mach numbers given as below. 


M ± = 


M±(M), 


(Wo) 
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if \M\ > 1, 

otherwise, 



where 


(146) 


Mf(M) = ±i (M ± l) 2 ± /?(M 2 - l) 2 , 

satisfy the properties [M1]-[M6]. .M^Lq reduces to the Van Leer formula [6] and is the 
lowest degree polynomial continuous at M — il. Figure 1 displays the distribution of M 
with various values of /?. 



- 1.0 0.0 1.0 

M 


Fig. 1 Mf vs M ; solid line: 0 = 0, ■ ■ • : 0 — —1/16, : 0 — 1/2, oo : 0 — 1/8. 


Let 


Mf +1/2 = \(M j+1/2 ± \Mj+i/ 2 1) = \M j+ i / 2 (1 + sign(M i+1/2 )), § (15) 

combining (6), (7), and (15) gives 

F S+V2 = ■ {16) 

Note that as I s identically equal to zero, the sign of Mj+i / 2 is immaterial since the 

convective flux vanishes with iWj+x/2* This also means the switching is continuous, but 
not differentiable. 

We remind the reader that the Mach numbers appearing in (13)-(14) need not be 
based on the speed of sound at grid dj,j € Z. In fact we propose in this paper to use 
a general interface value a i+ 1 / 2 , which will be determined in §2.3 to achieve some unique 
property. 

§ We caution the reader that here the “ ± ” signs in the interface Mach number have 
different meanings than that attached to the cell- point definitions in (14a,b). 
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We turn now to the pressure splitting in which the interface pressure is defined by 


Pj+ 1/2 = P~j Pj + Pj+lPj+ 1 * ( 17 ) 

Noting that we propose the use of normalized split pressure polynomials p ± : M E 
U —y p ± (M){ 0, 1]. First observed in reference [22], the pressure can be expanded using the 
same set of basis functions, (M ± 1), as clearly implied from the characteristic equations 
(9). 

Proposition 2.2: We require that the split pressures p ± satisfy the following properties: 
[PI] p + + p~ = 1, for consistency. 

[P2] 0 < p± < 1, as required by the physical constraint that the pressure be nonneg- 
ative. 

[P3] P + is a monotone increasing function of M, but p monotone decreasing with 
M, thereby allowing proper transition for upwinding. 

[P4] p + (M) = p - (— M) 

[P5] p + = 1 as M > 1; p~ = 1 as M < —1. 

[P6] p± are continuously differentiable. 

Similar assertion can be made about the pressure splitting as for the Mach number. 
We now prove the following properties for Pj+i/ 2 - 

Lemma 2.2: Pj+ 1/2 € [0,py +p>+i]. 

Proof. The proof follows directly by virtue of [P2] and 0 < |(1 i sign(M)) < 1. 

Remark: One can also show that pj+ 1/2 I s a monotone function of both p } and Pj+i, and 
monotone increasing in but decreasing in Mj+ 1 due to requirement [P3]. 

Definition of P: The split pressures given below, 

p±= Ji(l±s ig n(M)), (lga) 

[ pJ(M), otherwise. 

where 

Pa(M) = i(Jlf ± 1) 2 (2 t M) ± aM(M 2 - l) 2 , (186) 

hold [P1]-[P6]. Clearly, corresponding to Van Leer’s formula, are the lowest degree 
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differentiable polynomials satisfying [P1]-[P6]. The p* vs M curves are shown in Fig. 2. 



Fig. 2 vs M\ solid line: a = 0, • • • : a = —3/4, : a = —1/4, oo : a — 3/16. 


Next we suggest the values of parameters (a, 0) to be used. 

Lemma 2.3: To satisfy properties [M1]-[M6] and [P1]-[P6], we have 

-3/4 <a < 3/16, (19a) 

-1/16 </3 < 1/2. (196) 


Proof. The proof is simple and omitted here. 

Lemma 2.4: The parameters (a, /3) selected on the basis of the criteria: 


and 

d 2 M± 

n£w- 0 ’ 

<Pp °(±l) = 0 
dM 2 1 ' ’ 

are 

a = 3/16 and 0 = 1/8. 


(20a) 

(206) 

( 21 ) 


Proof. The proof is trivial and omitted here. 

The curves corresponding to these values are included in Figs. 1 and 2 respectively. 
These values are recommended for use as this pair have given the best performance over 
many problems tested [21], including those shown in the present paper. 
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Lemma 2.5: Given a bounded set K of U, then for every states (Ml, Mr) and (Ml, Mr) 
in K, there exists constants C ± (/?) > 0 and C ± (a) > 0, such that 


Similarly, 


\Mt-Mt\<C + mM L -M L \, (22a) 

IMg - MkI < C~(0)IMr - MrI- (22b) 

\pt-Pt\<C + (a)\M L -M L \ (23 a) 

\p- R ~p- R \<C-(a)\M R -M R \. (23 b) 


Proof. Referring to Fig. 1, we make use of the differentiability of M ± in (14). Setting 

C + (fi) = max M eTi\~^U 

_ . dM~ . 

C (p) = maxMeTC | -ffjyT I 

completes the proof. The proof for (23a, b) follows the same steps. 

Remark: Under (19b) in Lemma 2.3, namely —1/16 < (3 < 1/2, one easily gets C ± = 1. 
Lemma 2.6: The interface Mach number defined in (8a, b), 

Mlr = tti(Ml,Mr) — M+ + M^, 

is a Lipschitz continuous function of its arguments and satisfies the consistency condition 

m(M, M) = M. 

It is monotone nondecreasing in its first argument, nonincreasing in its second. 

Proof. The proof for Lipschitz continuity is equivalent to showing for every states (Ml , Mr ) 
and (Ml, Mr) in K of U that 

| m(M L ,M R ) - m(M L ,M R )\ = |M+ + M* - M+ - Mr\ 

= \(M+ - M+) + (M- - M~ r )\ 

< \(M l -M l \ + \M r -Mr\. 

The proof of consistency and monotone property is straightforward and omitted. 

Remark: One can assume in general for higher-order formulas that 

M j+ !/ 2 = m(Mj- k+ 1 ,...,M j+k ), 
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and similarly 


Pj+ 1/2 — p(Afj— jfc+lj •••■>Mj+k',Pj-k+h ’—iPj+k), 
for k > 1 and show the Lipschitz continuity and consistency condition. 

2.3 Definition of aj+ 1/2 

Recall that we still have one yet unattended task: namely selecting an interface speed 
of sound a,j+i/ 2 in order to achieve the unified formula for both the U- and M-splittings. 
To achieve this unification, it is obvious we can no longer use each respective speed of 
sound, a,j or o,j+ 1 , but instead should use a common one: 

«i + i/2=a(U;,U, + i). (24) 

In other words, the interface Mach number should be defined on the basis of a properly 
defined speed of sound. Since the interface flux is viewed as a recipient of contributions 
from both neighboring cells, it makes sense to use a common speed of sound evaluated 
there, as a basis for determining upwinding. It is this freedom that allows us to make a 
favorable choice to attain an exact capture of a single stationary shock. This concept of 
using a common speed of sound is also employed in [11]. 

Lemma 2.6: Consider a stationary shock problem: 

f Ul, x < Xj, 

\ Uk, x > xj , 

where the normal shock relation is connected by U l «md U/j and Uj > fly- Let fly be the 
critical speed of sound evaluated at Uy(= Ul), then 

Vj+i/2 = af/uj ( 25 ) 

exactly captures the shock without any intermediate states. 

Proof. Referring to Fig. 3, we consider at most one intermediate (numerical) state “ 1 ” 
between the “L” and “R” states. ( 


• 

m 

• 

L 

i 

R 


- 1/2 1/2 


Fig. 3 One intermediate shock cell between “T” and “R”. 
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Let us now express the numerical flux at interfaces (denoted here respectively by 
=pl /2) enclosing the intermediate cell i. In order to explicitly include the undefined speed 
of sound, we opt for using the velocity (U-) splitting. Similar to (8b) and (14a, b), (or 
referring to [2],) we write 

U-i/2 = ut + u 7 ■ 

Assuming Ml > 1 and M< < 1 without loss of generality, (since there must be a subsonic 
point connecting the supersonic point in the case of normal shock) then we have 

1/2 = U L — («» - a_ii 2 f g + {ui,a-i! 2 ), 

where we define 

g ± (u, a) = — (1 + 4£(^ ± l) 2 ). 

It is easy to show u_!/ 2 > 0 as ul > hence by (16) we get 


j m -l/2 \ 

F_i / 2 = u L m_ 1 / 2 +P-1/2 I, rn-1/2 = Plu-i/2- 

\ h t L m - 1/2 / 

This flux must be balanced with F l as the U L" state is fixed. Hence, a simple algebra 
yields 

( u -i/ 2 ^ = 

V P— 1/2 / \PlJ 

These two conditions can be satisfied exactly, for any (p»,Pi), by requiring only one condi- 
tion: 

Ui = a_ i/ 2 . (26) 

Let us now turn to the “1/2” face. The interface velocity becomes 

Ul /2 = (ttj + a 1 / 2 ) 2 ^ - (u,,a 1 / 2 ) — (ur - ai/ 2 ) 2 g + (uR,a 1 / 2 ) > 0. 

The inequality holds as Ui > ur is expected. Again using (16) yields 


m 1/2 

F 1/2 = [ u i m 1 / 2 +p 1 / 2 J , 
ht t m 1/2 


”*1/2 = 1/2- 


In order to eliminate this intermediate state, we let 


U,- = Ur, 
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It turns out that the last two requirements in the above equation respectively for p t and 
Pi are automatically satisfied if (26) and the first requirement, ur = u,-, are set. Thus, 


0—1/2 =UR = Ui . 


This requirement in fact enforces the conservation of fluxes by relating the upstream super- 
sonic U L" state to the downstream subsonic “R" state. It is desirable to express a_i /2 in 
terms of the upstream state “T”. In fact, the well-known Prandtl relation (see for example 
[23]) is at hand for use. 


uru l = a * L 2 


_ * 2 _ 
— a R ~ 


7+1 


a\ 


where a* is the speed of sound based on the total enthalpy hf . Putting L — > j and R * j '+1 , 
we complete the proof. 

Remark: It is well known that insofar as determining whether an isoenergetic flow lies at 
the sonic or in the supersonic or subsonic regime, the critical speed of sound defined above 
can be used as a reference [23] . 

Note that the formula (25) is valid for ul > a* L . For a flow covering different speed 
regimes, we extend the definition to: 


aj+ 1/2 — min(a£, a#), 


a 


a = a*min(l, — ). 

\u\ 


(27) 


That is, a is taken to be a* when |u| < a*. 

Remark: Other formulas can be also used for simplicity, but at the expense of losing the 
above exact property. We have tested the following obvious alternatives: 

a j + 1/2 = 2 a j + 1 ) ’ 

a j+ 1/2 = y/ a j a j+l i 
Oj_|_i/2 = mm(oj , Uj+i). 

Despite some expected minor differences among their solutions, no other adverse effects 
have been discovered thus far in our tests. 

Remark: The above exact shock capturing property holds for any (a, /?). 

Remark: The Van Leer splitting modified by the above procedure also allows exact reso- 
lution of a shock, see later in Fig. 5. 

Remark: Due to the above exact property, the AUSM + , like the Roe splitting, also can 
not distinguish whether the discontinuity is a shock or an expansion. However, it must be 
noted that except in this rather isolated situation, the AUSM + does not behave like the 
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Roe scheme in other cases where the latter is known to yield unsatisfactory solutions. (See 
examples 6 and 7 in section 7.) More specifically, the AUSM + , based on our numerical 
experiences, does not appear to branch into an entropy-violating solution if the initial 
condition is not exactly the expansion-jump condition. 

In s umm ary, the above splittings for both the advective and pressure terms completely 
define the Euler numerical flux. The AUSM + algorithm can be simply summarized as 
follows: 


Let j and j + 1 states be given, then 

(1) : Mj = Ujjo,jj r \j 2 i and Mj — u j+i j Oji-i-i j 2 via (27), 

(2) : Mj+i/ 2 = Mj + M ~ +1 , and Pj+ 1/2 = pjpj +Pj+iPj+i, 

* / 


(3): Get Mf +l/2 via (15) and F j+1 / 2 = a j+ 1/2 < 

+ Pj+1/2 





+ M j+ 1/2 



j + 1 / 


3. ANALYSIS OF AUSM+ 

Definition 3.1: The numerical flux defined above is rewritten as, 

F y+ 1/2 = Mj+i/ 2 -(&j + &j+ 1 ) - -\Mj + 1 / 2 \A j+1 / 2 $ + Pj+ 1 / 2 , ( 28 ) 

where (M, p) J+ i/ 2 are defined in sections 2.2 and 2.3, and Aj + i/ 2 {*} = — {•}>• 

While the associated properties of the interface quantities (M,p, a) J+1 / 2 have been 
discussed in the last section, we will address in what follows the mathematical consequences 
of the complete flux formula (28). 

Remark: The first term on the RHS is clearly not a simple average of “j” and “j + 1” 
states, but rather a Mach number-weighted average. 

Remark: The dissipation coefficient \M j+1 / 2 \ is merely a scalar — hence the system is 
decoupled after M, +1 / 2 has been defined, thus requiring only 0{n ) operations, in contrast 
to 0(n 2 ) operations by the other FDS. 

Remark: The present scheme does not involve differentiation, specifically the Jacobian ma- 
trix, in the evaluation of F J+1 / 2 ; it always involves only the common term M j+1 / 2 for any 
additional conservation laws. Hence, it is readily extendable to a general equation of state 
and non-equilibri um flows. Again, the cost is only linearly increased with the additional 
conservation equations considered. In the case of the Roe or Osher scheme in which the 
Jacobian matrix is essential, there is a need for redefining the averaged/intermediate states 
as additional equations are included. 
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We now proceed to prove some interesting numerical properties of the present scheme. 
For the sake of algebraic clarity, we choose to use variable u, instead of M, in the analysis. 
Recall that they are interchangeable, for j € Z, 

u j+ 1/2 = i/ 2 a i+i/ 2 - 

Lemma 3.1: The present sphtting preserves stationary contact discontinuity. 

Proof. Let pj = p j+ 1 = p, Uj = itj+i = 0, and pj / p J+ i across the stationary contact. 
Then M j+1/2 = 0, Pj+ 1/2 = P, and hence F^ 1/2 = 0. That is, no non- vanishing numerical 
convective flux is created at the interface. Hence, the flow remains stationary and the 
contact discontinuous. 

Lemma 3.2: For steady flows, the present scheme preserves the constancy of total en- 
thalpy. 

Proof. Assuming without loss of generality a unidirectional flow, e.g., u > 0, the discrete 
continuity equation yields for steady flow, 

u j+i/2^i — u ^-i/2Pj~ 1 = 

And the energy equation gives 

u "j+i/2 Pfitj ~ u j—i = 

Hence, we find htj = htj - 1 = • • • = hto, and the proof is complete. 

We note that this property is not satisfied by many upwind schemes, including those 
by Godunov, Roe, Osher, Steger- Warming, and Van Leer. To my knowledge, only those 
which upwind the unsplit form of the total enthalpy will preserve this quantity [1,2,11,15- 

17]- 

Lemma 3.3: Under the CFL-like condition, 

0 < A «+ +1/2 + (-u“_ 1/2 ) <1, A = At/ Ax, (29) 

then (a) the density function is monotone, and (b) the scheme preserves positivity of 
density. 

Proof. The mass conservation yields for j € Z, 

p’ +1 = p? - x (Kh- i/ 2 
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where 

fp,. 7+1/2 = “7+1/2 Pj + U j+l/2pj+l- 

Substituting and rewriting, 

p] +1 — ^“7-l/2^7-l + 1 _ ^(“7+1/2 - “>- 1 / 2 ) Pj + (~^ U j+l/ 2 ^Pj+l 

= ( 30& ) 

Since u + > 0, and u~ < 0, / is a monotone function of its arguments Pj+i) iff 

(29) is true. 

Next since > 0,Vj, (30b) immediately gives 

/>" > 0, n > 0, Vj, 

under the condition (29). Hence, the positivity of density is preserved. And (29) is called 
the positivity condition. The issue of preserving positivity of species mass fractions has 
also been studied by Larrouturou [24]. 

Remark : For linear case, Uj— 1/2 = Uj+ 1/2 = “ =const., it is easy to show that (29) 
satisfies the TVD condition of Harten [25]. But for nonlinear case, condition (29) does not 
guarantee TVD. 

Remark: Let 

u t - M t a 3 + 1/2’ “7 - a j-i/ 2 i 3 € Z , (31a) 

then we have the following inequality, 

<+l/2 + (- U 7-l 12) < u t~ U J- ( 316 ) 

Proof. By definition, 

2u 7+1/2 = U i+ !/ 2 + K+I/ 2 I = “7 + “7+1 + l“7 + “7+1 1 
2 u 7-!/2 = U j~ 1 / 2 ~ K-l/zl = u t - 1 + “7 ~~ l“7~l + u j I 

Then 

2(“t+ 1/2 - “7_ 1/2 ) = “7 + “7+i - u U ~ “7 + 1 “7 + u 7+ii + l“7-i + “71 
< “7 + u 7+i ~ “7-1 - “7 + u 7 _ “7+i + “7-i _ “7 

= 2(«7-«7). 
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Lemma. 3.4t The above CFL-like positivity condition can be expressed in terms of a 
stronger condition for a given A, 

0 < A uj + (—uj) < 1. (32) 

Proof. Using (31), it is easy to see that (32) guarantees (29) for a given A. 

Comparison of Well Known Numerical Fluxes 

Next, we write well-known numerical fluxes in a form similar to (28) for comparison. 

Roe: 

F i+1/2 = 1 (w*)j + («*),•+ 1 - 1 |A(U) 1/2 | Aj +1 / 2 U + -(Pj + Pj+i), (33a) 

where U is the Roe- averaged state, and ]A| = A~*~ — A in the usual sense. 

Osher; 

F y+1/2 = \ [(«#) f + («*)/«] - \ |A(U)| dU + i(P i + P>+l). (336) 

Steger- Warming: 

Fj+i /2 = 1 + (u $) J+ 1 - 1 A i+1 / 2 |A|U + -(Pi + P;+i)- (33c) 

Van Leer/Hanel: 

Fj +1/ 2 = \ (M*)j + (M*)y+ 1 - 1 A i+1/2 |M|$ + Pj+i/2* (33d) 

Remark: The Van Leer/Hanel splitting is seen also to require only 0(n ) operations, since 
\M\ is a scalar. 

Remark: The Steger- Wanning and Roe splittings differ only in the evaluation of the ab- 
solute Jacobian. The absolute quantity in the latter is evaluated using both the “j” and 
u j -}- 1” states and is taken outside of the difference operator. In fact, this comparison 
reveals clearly that a striking difference in form between the FVS and FDS lies in whether 
the dissipation matrix ( or scalar) is differenced. 

Remark: In the above sense, the present scheme may appear formally close to a FDS, but 
it differs in the averaged term. On the other hand, the method retains the efficiency of the 
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Van Leer scheme in defining the dissipation term. Moreover, it achieves the high level of 
accuracy attributed only to the Roe and Osher methods, as will be borne out by a variety 
of Euler and Navier-Stokes calculations. Consequently, the present scheme is neither FDS 
nor FVS, but rather a hybrid one. 

Remark: It is of interest investigating the limiting form of the flux valid for the boundary 
layer flow where the transverse Mach number is small, M x / 2 « 1. Assuming M x j 2 > 0, 
the y-momentum flux (cf. (28)) becomes, by retaining only leading terms: 

G j+1/ 2 = M#(l-^A i+1/2 M)-^M(l-iA i+1/2 M)A$+P(l-^A i+1/2 pAf+---), (34) 

where / = \(fj + fj+i)- Since the leading terms coincide with the central-difference 
approximation, the nominally first-order upwind formula (34) now tends to be second- 
order accurate as M — * 0. 

Variants in Numerical Energy Flux 

Furthermore, we propose another viable variation which we have found to be reli- 
able. In the energy equation, the term (up) can be separated from the total flux puh t . 
Specifically, we suggest the following two options. 

(puh t )j+ 1 / 2 = u ~j +x / 2 Pj e tj + U J+ 1/2 Pi+ 1 e< i-i- 1 + u i+i/ 2 Pj+i/ 2 - (35a) 

The terms u j+x / 2 and p j+ x / 2 are defined using (8b), (17), and (27). 

Or following [22], we can derive 

(puht)j+\/ 2 = u ^+i/ 2 Pjhtj H" u j+i/ 2 Pj+ihtj+i + 4<7a J+ i/ 2 (Af^ Mj — Mj +l Mj +x )., (35 b) 

where the parameter a gives reise a family of choices. As a = 0, we recover the basic total- 
enthalpy-preserving scheme recommended in this paper. The original Van Leer splitting 
sets a = 2h+a 5 ' Note that the last parenthes in (35b) can take either positive or negative 
value although M~ < 0. 

4. HIGHER-ORDER EXTENSION 

While the first-order schemes are simple and useful for analysis, it is generally desirable 
to employ higher-order schemes. The extension, although aiming at improving accuracy, 
should as well preserve those properties arrived at for the first-order scheme, such as 
monotonicity, positivity, efficiency, and robustness. We will attempt to focus on these 
properties in this section. 

The MUSCL [26] strategy for preserving monotonicity is adopted here, via the use of a 
nonlinear limiter sensing the ratio of neighboring first-differences of appropriate variables. 
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Many forms have been proposed in the literature and to our experiences they all work 
reasonably well although certain differences in computed results can exist. Systematic 
discussion on the subject of limiters is beyond the scope of this paper. We onlj give 
the formulas used in our calculations, acknowledging that there may be a more elaborate 
choice. In [21], we will address these issues in more detail. 

Following the spirit of the present method that the fluxes can be evaluated sequentially, 
we avoid using the characteristic variables in the extrapolation procedure. It may suffice to 
extrapolate the primitive variables W, thereby leading to a simple and efficient procedure: 

W[ := W, + l^r,)(Wy - Wj-j), W = (p,u,pf, 

W« := W i+1 - ltf(— )(W j+2 - Wy +1 ), (36) 

2 r j+ 1 

where the subscripts u L ri and represent the states that take the place of j and j -1-1 
in the formulas derived in section 2 for the first-order scheme. To take grid nonuniformity 
into account, grid-size weighting can also be formulated into (36) if desired, see for example 
[27]. The function ip is a limiter whose argument is defined as 


_ A i+ i/2^ 

3 A j-i/ 2 Wi ' 


W = (w u w 2 ,w z ) T = (p,u,p) T . 


(37) 


The minmod and superbee functions are normally chosen in our calculations. The effects 
of limi ters on either accuracy or convergence rate axe often present, but discussion of them 
is beyond the scope of this paper. A note of suggestion may be warranted in the case 
of viscous calculations. For the reason of maintaining uniform accuracy and convergence 
rate, it is desirable that the limiters be set to unity insofar as stability is allowed, so that 
physically true extrema (such as the peak temperature in a thermal boundary) will not be 
clipped or the limit-cycle behavior can be minimized. 

In what follows, we will show that the second-order scheme satisfies monotonicity and 
positivity conditions. 

Let us again consider the continuity equation written in the finite-volume form (30a), 


Pj +l = Pi 


~ A (fp,}+ 1/2 fpj~ I/ 2 ) ‘ 


(30a) 


Now we write the mass flux at u j + 1 /2” as 

fp,j+ 1/2 = U ~j+ 1 / 2 PL u j+i/ 2 PR‘> 


(38) 


where, by virtue of (36), 

Pl = P] + /2 ~ Pi- 1) = P L (Pj-i ’ Pi ’ ^i+1/2 )’ 

Pr = P ]+ 1 — 2 ^;+ 1 / 2 - Pi+ 1 ) = P* 0 >] + i,PH 2 , 1 > i+ x/ 2 )- 
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(39a) 

(396) 


In the above expression, we have adopted new notation for the limiters, for algebraic clarity. 


^ + +i/2 = ^( r ;)’ ^i+i/2 = ^(^7)- ( 40 ) 

Note that the superscripts “+” and are used to associate them with forward and 
backward extrapolations, but without attaching sign values to them. Using (36), we express 
the limiters as 

^y+l/2 = ^j+1/2 ~ ^j+l/2^Pj iPj+l'Pj+2)’ (41) 

and hence the interface densities as 

Pl ~ PL(Pj-ii p ] ? Pj+i)i p'r = PRiPj ’ Pj+i* Pj+2)‘ (42) 

Similarly, we can write 

u j+i/2 — u j+i/2( u y-n u y i M i+i» u i+2)- (43) 

That is, we have a 5-point scheme in (30) as expected for the second-order method. By 
substitution, 

p" +1 = p"-a{u + +1/2 [(.7 + 5^ +1/2 (/>7 - />"-!)] + “ 7 + 1/2 [/>”+. - 5 ^ 7 +l / 2 «+2 - P?+l) 

- “7-1/2 ^i-l ■*" ^i-l/z^i-l ~ Pl- 2 ) “ “i- 1/2 Pj ~ ~ }' 

(44) 

To utilize the simplicity found in the proof of Lemma 3.3 for the first-order scheme, we 
first group the basic first-order terms, and then rearrange the remaining terms by utilizing 
(37), so that the resulting equation appears formally like the first-order scheme as in (30b). 

After some algebraic manipulations, we get 

P] +1 ~ P n j ~ ^ j u >+i /2 P] + u j+i/2P 1 j+i ~ U /-l/2^i-l “ U j-1/2P] | 

+ 2 | - “i+i/2^i+i/2(^i “ Pj-i) + u j+i /2^j+i/2 r i+ 1 (Pj+i ~~ Pi ) 

+ U ^-l/2 ~~~^~(Pj ~ ^i _1 ) — U j-l/2 X ^j-l/2 r j+ 1 (Pj+l ~ Pj )|- 
Or 

p” +1 = Ci-ip?-! + C jP ] + c j+1 p] +l , (45) 
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where 


Ci = 


X 

X "2 


Cj-i 

c j+ 1 


«; +1/2 ( 2 + ^ + +1/2 ) + (-«-_ 1/2 )(2 + V’-_ 1/2 ) 

/- + ^-1/211 
+ U i+l/2 r i+l^i+l/2 _u i— 1/2 rj _ 1 Jj’ 

2 | U i+l/2^/+l/2 + U y’-l/2( 2 r> _ 1 )}» 

2 U J-l/2)^y-l/2 + U i+l/2)( 2 ~ r J+lV’_ ?+ i/ 2 )|- 


(46a) 


(466) 

(46c) 


Definition 4.1: The second-order accurate density function in (45) and (46) is monotone 
\/j and satifies the positivity condition Vn > 0 iff 

Cj -1 > 0, Cj > 0, and, C j+1 >0, j € Z. (47) 


Let us define 


and 


u min < u j-l/2i u j+ 1/2 ^ Wmai) Vi, 

u mar = maX [“max, 0] , U~ in = min [u m .n, 0] . 


(48) 

(49) 


We turn now to show the conditions for the limiter function and a CFL-like inequality. 
Lemma 4.1: Under the following conditions, 

0 < r/j(r) <2 r < V’max, (50) 


and 

- »^.X1 + ~) < 1. (SI) 

then (a) the density function is monotone, and (b) the scheme preserves positivity of 
density. 

Proof. Since ut ±1/2 > 0, and uJ ±1/2 < 0, and requiring 


$j± 1/2 — 0 » 


we have the following inequality, 


Cj-i > 


ut 


i- 1/2 /n ^i- 1/2 


r,-i 


-( 2 - 


)> 


( 52 ) 
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and 


'i + 1 ^ 


“i+l/2 


(2-r J+1 0 i+1/2 ). 
This immediately yields that Definition 4.1 (47) is met iff 

1/2 


r i-i 


r i-i 


< 2 , 


r }+ l^y+j/2 = r >+ 1 ^^') - 2 ’ 


and 




U i"-fl/ 2 (^ + ^j+ 1/2) + ( U >-l/ 2 )(^ + ^'-1/2) 




+ U i+l/2 r >+l^j+l/2 U T-l/2 3 r .J x ] } 
Using (48), (49), and (52), one can easily show 


> 0. 


Cj > 1 


U, 


:(2 + ^1/2) + ( _u min )( 2 + $ j- 1/ 2 ) 


(53a) 

(536) 


(53c) 


Assuming 0 < ifrf ±1 / 2 - Vmax, the above equation leads to 

K U max ~ U min)( l + ^—) < l - 

This completes the proof of (a). The proof of (b) follows immediately as in Lemma 3.3. 

Hence, we see that the positivity of density can be preserved from the “second-order 
(except at extrema)” scheme, (30a), (38), and (39) under conditions (50) and (51). 

Remark: Similar proof can be found also in the TVD-type study. Also a weaker condition 
than (51) may be possible, e.g., by admitting < 0, hence maintaining second-order 
accuracy at extrema. However, this is not the purpose of the above proof and is beyond 
the scope of the present paper. 


5. GENERALIZED CURVILINEAR COORDINATES 


This section gives a prescription for treating generalized curvilinear coordinates. To 
facilitate the description, let us first define the notation for the relevant variables in the 
3D Euler equations. The physical variables in a phase space of dimension 5 are denoted 
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by a boldface uppercase letter or column vector whose elements are denoted by lowercase 
letters. 


U = 


p \ 


f P \ 

pu 


pu 

pv 

, * = 

pv 

pw 


pw 

\pe t / 


\ph t / 


(54) 


where e t = e + 0.5(u 2 + v 2 + w 2 ) and h t = e t + p/p ■ The geometrical vectors in physical 
(Cartesian) space of dimension 3 are denoted by an overhead arrow “ . The fluid velocity 


is 


V = ui + vj + wk, 

and the normal vector of the boundary surface of a control volume 

S — SjZ “I - SyJ “f S z k. 

The inviscid fluxes in 3D physical space are compactly written as 


= V + P = F^ + P, where P = 


(55a) 


(555) 



P \ 


/°A 

F = 

pu 


pi 

pv 

PI 


pw 


pk 


\phtJ 


\0/ 

The first term in F 

is the flux of 


/°A 

pi 

pi 

pk 

\0/ 


(56) 


simply the pressure flux. 

The inviscid conservation laws can be conveniently expressed over an arbitrary control 
volume Q, in an integral form: 


/ 


^-dv+ f[&V + P].dS=0. 
ut J 
9. an 


(57) 


The discrete form, describing the rate of change of U in Vij t k via balance of fluxes 
through all enclosing faces, 5/, l = 1, • • ■ , LX(i, j, fc), can be cast as 


A , LX 

u "S = u "m-— E(f! c) +P')*^=°- 

j ^ ^ 


(58) 


In what follows, we give the algorithm for constructing the numerical fluxes in an 
arbitrary finite volume, based on that described in Section 2. Let us consider the situation 
in which the computation domain is divided into hexahedrons where each pair of opposing 
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faces constitute the coordinates £,77, and ( with corresponding discrete indices i,j, and 
k, respectively. For the purpose of presenting the procedures, it suffices to consider one 
interface in the ^-direction, whose surface vector is being 

the unit normal vector and S l ^i j 2 ,j.k the area. Again, we find it convenient to use the 
convention of the “L" and “R” states for representing conditions on each side of the 
interface Sf +1/2J k . 

(1) Project the velocity vector at the cell centers (i,j,k),(i + 1, j, k) to Sf +l ^ 2 - k , 

Vl = Vi,j,k • n$, Vr = Vi+i'i^k • n$. (59) 


Hereinafter we use “overhead” ( ) to denote quantities in curvilinear coordinates. 
It is reminded that if higher order accuracy is desired, extrapolation given in 
section 4 can be adopted for V{j t k and Vi+ 

(2) Define the corresponding Mach numbers, 


M l = 


Vl 

<H+\/2,j,k 


M R = 


Vr 

ai+i/2,j,k 


with a suitable definition of Ui+ i/2,j,k such as (27). 

(3) Define the interface convective Mach number by writing 


M i+ l/2,j,k — Mr + Mr , 


(60) 


(61) 


where 


m£ = 1 2 


\{Mi + \Ml\), if \Ml\ > 1, 


( 1 


Mt ( M l ) , otherwise, 


Mr = { 


-(Mr-\Mr\), if \Mr\ > 1, 
MJ (Mr), otherwise. 


The formulas for Mf(M) are given in (14b). 

(4) Obtain the interface pressure Pi+i/ 2 ,j,k in a similar fashion by using the above- 
defined Mach numbers in (17)-(18), 


P»+l/2,i,fc — Pi+l/2,j,k 


(l) 


5 : 

\ 0 / j +i/2J,k 


(62) 


(5) Assemble the interface numerical flux, 

(F • S) i+ i/2,j,k = {-^5.1/2,;,* + M- +1/2J ' k $i+i,j,k + Pi+l/2,j,fc} Si+l/2,j,k- 

(63) 
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Lemma 5.1: The CFL-like condition of the first-order scheme for the 3D case is 


0 — ^i-l/2,j,k^ 1 k/2,j,k)~^ 

u i,j,k 

(^Cj,k+l/2^i,j,k+l/2 ~ - 1? 

as the density satisfies the monotone condition and remains positive. 

Proof. The proof directly follows that in Lemma 3.3. 

Remark: A similar formula for the second-order accurate scheme also can be derived, see 

Lemma 4.1. 


6. BOUNDARY CONDITIONS 

This topic generally does not get enough space in a paper. From an analytical point of 
view, the boundary conditions required to close the mathematical system are well defined. 
When these analytical boundary conditions are transformed in the numerical discretization, 
ambiguities and difficulties begin to compound with the interior flux schemes. While 
a number of papers have been devoted to proposing and analyzing numerical boundary 
conditions, unfortunately this issue still has not received sufficient investigation. Often it 
is barely mentioned in a paper, and the exact implementation, where variations begin, is 
generally left out. The situation is even more disappointing in the case of unsteady flows. 

We also assert that there is a tight relation between the interior scheme and the 
selection of apppropriate boundary conditions. In other words, the concern of compatibility 
between a numerical flux scheme and boundary scheme must be observed - a numerical 
boundary procedure working properly for scheme A may not be suitable for scheme B. To 
our knowledge, little has been reported on this aspect. 

Of various types of boundary conditions, the treatment at a physical wall (solid or 
porous) perhaps is one of the most diverse among CFD practitioners, and is the one 
that affects solution behavior (accuracy and convergence rate) the most if due care is 
not exercised. Although this situation has been observed, to sort out the issues involved is 
beyond the scope of the present paper. Nevertheless, we give in the following the procedures 
used for inviscid solid wall conditions, on which all comparisons among methods are based. 

A linear extrapolation of all variables from the nearest interior points (two points 
in the case of second-order solution) to the wall is employed. For slip wall, we impose 
the condition of zero velocity component normal to the wall to recover the Cartesian 
components. In Fig. 4, we let 


Vet = Ul + J-(Vl - V 2 ), 
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and 


V n — ^ei ' Tlw- 


Then a simple vector relation yields 

V w — V ex V n n w . 

As a result, the inviscid numerical flux is simply 


F XV F yj • Tlw 


/ 0 \ 

Pw'R'Wx 
Pw^wy 
PwKwz 

\ o / 


where p w is extrapolated from interior point(s). 



Fig. 4 Wall boundary condition. 


7. RESULTS AND DISCUSSION 


( 65 ) 


( 66 ) 


In this section we demonstrate the capability of the proposed numerical flux by con- 
firming the mathematical consequences established in previous sections for various relevant 
problems. The problems are listed below and will be discussed accordingly. 

1. Stationary shock discontinuity, 

2. Stationary contact discontinuity, 

3. Hypersonic conical boundary layer, 

4. ID shock tube problems, 
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5. 2D channel flow, 

6. Blunt body flow, 

7. Odd-even grid perturbation problem, 

8. Shock diffraction around a corner. 

For validation purposes, analytical solutions will be used as much as they are available. 
Otherwise, solution will be compared with that obtained by using Roe’s splitting as it has 
been established as an accurate and perhaps the most popular upwind scheme in use today. 

The first three problems are ID steady flows, dealing with the capability of accurately 
resolving shock, contact discontinuity, and a thin viscous layer; the last one of these tends 
in the limit of Re -* oo to behave as a contact discontinuity. In problems 4, we are 
concerned with accuracy in calculating unsteady flows. In problems 5 and 6, we consider 
2D steady inviscid flows. Finally, 2D unsteady solutions are discussed in problems 7 and 
8 . 

In some problems, first-order solution is included. Test on first-order solution is very 
meaningful for it reveals the sheer accuracy of a scheme, and the monotonicity. It is also 
fundamental because the TVD schemes and higher-order solutions can be built up from 
the first-order scheme. 

Figures 5 and 6, respectively, show the numerical result for a single stationary shock 
and contact discontinuity. The present AUSM + , like Roe splitting, is designed specifically 
to produce exact solution in these cases. In contrast, the AUSM and Van Leer splittings 
resolve the shock with two intermediate points. As mentioned earlier, a specific choice of 
speed of sound (Eq. (25) or (27)) used in the Van Leer splitting can also yield an exact 
shock solution. Unfortunately, this modification is not sufficient to rescue the Van Leer 
splitting from producing “fatty” dissipation in the limit of vanished flow speed. Figure 6 
shows that both Roe and AUSM+ (also AUSM) splittings give rise to exact solution of a 
stationary discontinuity, while the Van Leer splitting increasingly smears the discontinuity 
as calculation continues. 

The quasi- ID conical flow at Moo = 7.95 over a cone of 10° half angle was calculated 
using 64 cells. The profiles of pressure, temperature, and transverse velocity component 
on 64 cells are shown in Fig. 7 for the first-order results and in Fig. 8 the second-order 
results. The Roe splitting and the AUSM + give nearly identical results; interestingly the 
convergence rates also show remarkable resemblance. It is noted that the first-order result 
is already as good as the second-order one. To further raise curiosity, we include a first- 
order coarse-grid (32 cells) result of AUSM + for comparison; the accuracy is astonishing in 
that there are only about four (4) cells for resolving the extremely high-gradient boundary 
layer and the shock is nearly captured exactly. Roe solution behaves identically but is not 
included for the sake of clarity in plots. 

Next we consider the standard shock-tube problems. Here we used 100 cells and the 
minmod limiter for the “second-order” calculation. The CFL number was set 0.45. For the 
Sod problem, Fig. 9 clearly demonstrates the comparable accuracy of the present and the 
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Roe solution. The second-order solution shown in Fig. 10 again confirms the accuracy. We 
also add the more difficult problems studied by Yee [28], corresponding roughly to cases 
C and E in [28] in terms of density and pressure ratios and initial Mach number except 
without chemistry effects. The former involves a much stronger contact discontinuity than 
the Sod problem; the latter has a shock moving slowly against a high pressure region, 
with characteristics similar to that studied by Roberts [29], and with an additional large 
density jump. Figure 11 begins to show the differences between results obtained using the 
AUSM + and Roe splitting. Immediately behind the shock, the AUSM + appears to give 
better agreement with the exact solution of velocity. In Fig. 12, a long wave oscillation 
of the Roe solution, observed by Roberts [29] and Lin [30], also appears in this problem. 
Both AUSM + and AUSM results seem to be free from this behavior, but AUSM gives rise 
a glitch in the solution due to a large gradient at the contact discontinuity. 

We turn now to 2D calculations. First we show, in Fig. 13, the Mach contours of an 
M -- 1.8 flow over a 15° ramp. We used the minmod limiter and CFL=0.8 on 138 x 38 
cells. The ramp shock is sufficiently strong so that a Mach reflection appears on the top 
wall. A similar but weaker Mach stem also shows up on the bottom wall. In both cases, 
a distinct shear layer emanating from the triple point can be seen. Detailed distributions 
along the top and bottom walls are dispalyed in Figs. 14 and 15, among AUSM, AUSM + , 
and Roe solutions; they essentially agree with each other, results by AUSM and AUSM + 
nearly indistinguishable. The convergence history in Fig. 16 in terms of £2 norm exhibits 
that the Roe solution stalls after having dropped about three orders, while the other two 
residuals continue to decrease, also with nearly identical pace. We stress that the presented 
solutions have been confirmed to be grid- independent, hence the fine-grid solutions are not 
included. 

Second case in this category is a standard 4% bump in a channel with Moo = 1-4. 
The results were obtained using 132 x 68 cells, minmod limiter, and CFL=0.9. Again, 
the profiles on the bump wall in Fig. 17 show that three solutions agree well and all give 
sharp representation of three shocks, two respectively located at the leading and trailing 
edge of the bump and one located further downstream and produced by the reflection of 
the leading-edge shock from the other wall. The convergence history of the Roe solution 
behaves similarly to the previous case, as displayed in Fig. 18, except that it stalls at a 
much lower level. The AUSM and AUSM + , however, continue to behave well. 

Let us now consider a supersonic = 6.0 flow over a circular cylinder. This 

seemingly benign problem in fact fails the Osher and Roe methods [21,30]. The so-called 
“carbuncle” phenomenon by Roe’s method is exhibited in front of the cylinder, as compared 
with the AUSM + result in Fig. 19. Figure 20 displays the solutions by AUSM and AUSM + , 
showing a crisp shock resolution at the centerline. In particular, the AUSM + exhibits the 
evidence of improvement over the AUSM by eliminating the post-shock overshoot and 
almost getting rid of numerical shock point. Figure 21 compares the convergence history 
of various solutions, showing the effects of order of spatial accuracy, grid size, and CFL 
number. It is surprising to see that the second-order method in fact converges faster than 
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the first-order method, even with the shock. Finer grid makes the convergence slower, as 
consistent with other published results. 

Next, we study a very interesting and benign problem, first reported by Quirk [4]. 
A plane shock is moving into a quiescent region in a long constant-area channel. The 
computation grid is perturbed at the centerline by a very small magnitude, ±10“ 6 in our 
case, alternately at odd and even points. The Roe solution, seen in Fig. 22(a), quickly 
develops an odd-even decoupling behavior, eventually leading to an unacceptable solution. 
The present method, on the other hand, still preserves the plane shock even after a long 
time, as shown in Fig. 22(b). The variables profiles along the centerline, corresponding to 
the contours, axe displayed in Fig. 23 for comparison — the AUSM+ gives monotone and 
well behaved solution. 

The last problem is the diffraction of a supersonic moving shock over a 90-degree bend 
and the resulting complex flowfield. Quirk [4] has painstakingly shown the complexity of 
the flow using grid refinement to resolve the fine details. He pointed out that there was 
some n um erical difficulty encountered in the use of the Roe splitting. Thus it would be 
an interesting problem to test whether the present AUSM + would have the same difficulty 
or even face another. First we show the density contours of first-order AUSM , Roe, and 
Godunov solutions on a 71 x 71 grid, t which is extremely coarse for revealing any intricacy. 
Nevertheless, it confirms the findings of Woodward and Collela [31] and Quirk [4] that the 
Godunov scheme, as formally extended to 2D in the same manner as the other schemes 
reported in the present paper, can yield discontinuous expansion fans as seen in Fig. 24(a) 
— this, apparently violating entropy condition, clearly causes concern. Figures 24(b) and 
24(c) show that both the AUSM+ and Roe splittings respectively are able to break the 
expansion fans, while the AUSM+ solution appears to have slightly more diverged fans 
than the Roe solution. A fine grid (400 x 400) calculation was made using minmod limiter 
and CFL=0.4. The result is depicted in Fig. 25. The intricacy associated with this flow 
is astonishingly rich, although what is contained at this grid resolution is still not nearly 
comparable to that captured by Quirk [4]. 

8. CONCLUDING REMARKS 

In this paper we presented the construction and analysis of a new numerical flux 
scheme — AUSM + . We gave the associated mathematical properties and proved its mono- 
tonicity, positivity, and CFL-like condition. We proposed a modification to the Mach 
number and pressure polynomials used in the AUSM. Furthermore, a judicious choice of 
numerical speed of sound for the interface Mach number was shown to give exact resolution 
for an ID steady shock and significantly improved shock resolution in 2D cases. The reli- 
ability of the new scheme was evident as well in calculating some unsteady problems that 

f These results were obtained with the generosity of Dr. Yasuhiro Wada who kindly 
provided his code for producing not only solution but also plots. 
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have f Eii led pro mi nent flux schemes. We also stress that in addition to the demonstrated 
accuracy and reliability, the AUSM + requires computational effort far less than the flux 
difference splittings and only insignificantly more than Van Leer’s flux vector splitting. 
Moreover, since an exact linearization of the AUSM + flux can be derived, a consistent 
im plicit, scheme using these Jacobians is available. Results on this aspect will be reported 
in the near future. 
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Fig. 6 Stationary contact discontinuity. 












Fig. 11 Problem C; second-order solution. 






Fig. 12 Problem E: second-order solution. 


Fig. 13 Supersonic ramp problem; Mach number contours. 
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Fig. 14 Supersonic ramp problem; profiles on bottom wall. 
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Fig. 15 Supersonic ramp problem; profiles on top wall. 



Fig. 16 Convergence history for supersonic ramp problem; second-order solution. 
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Fig. 17 4% blimp problem; second-order solution on the bump. 



Fig. 18 Convergence history for 4% bump problem; second-order solution. 
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Fig. 19 Supersonic blunt body problem; Mach contours, left: Roe method, right: AUSM + . 
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Fig. 20 Supersonic blunt body problem; second-order solution. 
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Fig. 21 Convergence history for supersonic blunt body problem. 



Fig. 22 Odd-even gird perturbation problem; Mach contours, top: Roe method, bottom: AUSM . 
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Fig. 23 Odd-even gird perturbation problem; profiles along the centerline. 
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Fig. 24 Supersonic corner flow problem; density contours of first-order solution, top: Gudonov 
method, middle: Roe method, bottom: AUSM + . 
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Fig. 25 Supersonic corner flow problem; density contours of second-order (minmod limiter) 
AUSM + solution on 400 x 400 grid. 
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